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The development of implicit upwind algorithms for the solution of the three-dimensional, time-dependent Euler 
equations on unstructured tetrahedral meshes is described. The implicit temporal discretization involves either a two- 
sweep Gauss-Seidel relaxation procedure, a iwo-sweep Point-Jacobi relaxation procedure, ora single-sweep Point-Implicit 
procedure; the upwind spatial discretization is based on the llux -difference splitting of Roe. Detailed descriptions of the 
three implicit solution algorithms arc given, and calculations for the Boeing 747 transport configuration are presented 
to demonstrate the algorithms. Advantages and disadvantages of the implicit algorithms are discussed. A steady-state 
solution for the 747 configuration, obtained at transonic flow conditions using a mesh of over 100,000 cells, required 
less than one hour of CPU time on a Cray-2 computer, thus demonstrating the speed and robustness of the general 
capability. 


L INTRODUCTION 

In recent years, significant progress on develop- 
ing numerical algorithms for the solution of the gov- 
erning fluid flow equations on unstructured meshes has 
been made (refs. 1-7). This progress includes improve- 
ments in solution accuracy as well as computational effi- 
ciency. For example, upwind methods have been devel- 
oped for unstructured meshes which arc based on the local 
wave propagation characteristics of the flow, and con- 
sequently, produce highly accurate solutions (refs. 2,3). 
Most of these upwind methods, however, use explicit 
time-marching schemes to integrate the governing equa- 
tions in time to steady state. The explicit approach is 
computationally efficient when applied to meshes that are 
coarse, but the rate of convergence deteriorates signif- 
icantly when finer meshes are used. For cases where 
finer meshes are used, either a mulligrid strategy for con- 
vergence acceleration or an implicit temporal discretiza- 
tion which allows large time steps is required to obtain 
steady-state solutions in a computationally efficient man- 
ner. Implicit upwind solution algorithms for unstructured 
meshes in two dimensions have been reported by the au- 
thor in ref. 8. These algorithms arc similar to the point- 
implicit scheme of Tharcja et al. (ref 9), although the 


methods of ref. 8 arc fully implicit and not point im- 
plicit. The purpose of the paper is to describe the de- 
velopment of three implicit upwind algorithms for the 
solution of the three-dimensional time-dependent Euler 
equations. The implicit temporal discretization involves 
either a two-sweep Gauss-Seidel relaxation procedure, a 
two-sweep Point-Jacobi relaxation procedure, or a single 
sweep Point-Implicit procedure. The spatial discretization 
of the scheme is based on the upwind approach of Roe 
(ref. 10) referred to as flux-difference splitting (FDS). The 
FDS approach is naturally dissipative, and consequently 
captures shock waves and contact discontinuities sharply. 
Detailed descriptions of the three implicit solution algo- 
rithms are given, and calculations for the Boeing 747 
transport configuration are presented to demonstrate the 
efficiency of the algorithms. 

2 . EULER EQUATIONS 

In the present study the flow is assumed to be 
governed by the three-dimensional time-dependent Euler 
equations which may be written in integral form as 

Jt J QdV + 1 {Enz + Fny + Gn * )dS = 0 (1) 

n an 



where Q is the vector of conserved variables, and E*f\ 
and G are the convective fluxes. Equation (1) has been 
nondimensionalized by the freestream density and the 
freestream speed of sound. Also, the second integral 
in Eq. (1) is a boundary integral resulting from appli- 
cation of the divergence theorem, and n x , n y , and n z are 
Cartesian components of the unit normal to the boundary 
surface. 

3. SPATIAL DISCRETIZATION 

The spatial discretization is based on Roc’s flux- 
difference splitting which is herein implemented as a cell- 
centered scheme whereby the flow variables arc stored at 
the centroid of each tetrahedron and the control volume 
is simply the tetrahedron itself. Consequently, the spatial 
discretization involves a flux balance where the fluxes 
across the four faces of a given tetrahedron arc summed 
as 
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53 HAS = 53 ( E "r + Fn„ + On. )AS (2) 

m= 1 m=l 

where AS is the area of the face. The flux vector JI is 
approximated by 


upwind-biased interpolation for is determined simi- 
larly. Also the parameter k in Eq. (4) controls a fam- 
ily of difference schemes by appropriately weighting A_ 
and A+. On structured meshes, it is easy to show that 
k — — 1 yields a fully upwind scheme, k = 0 yields 
Fromm’s scheme, and k = 1 yields central differencing. 

In calculations involving upwind-biased schemes, os- 
cillations in the solution near shock waves are expected to 
occur. To eliminate these oscillations flux limiting is usu- 
ally required. The flux limiter modifies the upwind-biased 
interpolations for q~ and q+ such that, for example, 

<r = <lj + |[(i - ks)A_ + (1 + ks)A + ] (6) 

where $ is the flux limiter. In the present study, a con- 
tinuously differentiable flux limiter was employed which 
is defined by 


2A_ A+ T i 
A* + A* + e 


(7) 


where c is a very small number used to prevent division 
by zero in smooth regions of the flow. 


4. TEMPORAL DISCRETIZATIONS 


Ii = \\li(Q+) + H(Q~) - |/i|((? + -(?")] (3) 

where Q~ and arc the stale variables to the left and 
right of the cell face and A is the flux Jacobian matrix 
given by dll/dQ. Also the tilde and the absolute value 
sign indicate that the flux Jacobian is evaluated using the 
so-called Roe -averaged flow variables and the absolute 
value of the characteristic speeds. 

The left and right states, Q~ and are determined 
by upwind-biased interpolatioas of the primitive variables 
q. In three dimensions, for a given tetrahedron j, for 
example, the upwind-biased interpolation for q~ across 
the common face between tclrahedra j and k is defined 
by 

9” = m + ^[(1 - «)A- + (1 + k)A + ] (4) 

where 

A+ = q k ~ q } (5a) 

A. = qj - (5b) 

In Eqs. (4) and (5), q 3 and q k are the vectors of primi- 
tive variables at centroids j and h, respectively, and q lt 
the vector of primitive variables at node i (the node of 
tetrahedron j opposite to the face being considered), is de- 
termined by an invcrsc-distancc-wcighted average of the 
flow variables in the tetrahedra surrounding node i. The 


The temporal discretizations are of the implicit type 
and are generally derived by first linearizing the flux 
vector li according to 


y/ n+t =// n + (g) 

where dll/dQ is the flux Jacobian A , as discussed before, 
and A Q = (? n+1 - Q n . Linearizing both flux terms on 
the right-hand-side of Eq. (3) using Eq. (8), and ignoring 
the tilde on the flux Jacobian, results in 
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m = 1 

where 1 is the identity matrix, “vol” is the volume of the 
tetrahedron j , and A Q m is the change in flow variables 
in each of the four tetrahedra adjacent to tetrahedron j. 
Also in Eq. (9) and A~ are forward and backward 
flux Jacobians, respectively. For flux -difference splitting, 
the exact Jacobian A (derivative of the right-hand side of 
Eq. (3) with respect to Q) is too expensive to compute, 


A. S' 
(9) 
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and thus an approximate Jacobian is normally used. This 
is accomplished by constructing the Jacobians making 
use of the fact that the forward and backward Jacobians 
should have non-negative and non-positive eigenvalues 
(characteristic speeds), respectively. This is accomplished 
by expressing alternatively the Jacobians using similarity 
transformations such that 


yl + = RA+R ~ l A~ = RA~ R~ x (10) 


where A + and A are diagonal matrices whose diagonal 
elements are the eigenvalues A + and A” defined by 

= |(A + |A|) A~=i(A-|A|) (11) 

and R is the matrix whose columns are the corresponding 
eigenvectors. 

A similar implicit temporal discretization that is more 
efficient than the discretization of Eq. (9), is derived by 
linearizing the fiux vector of the quasi-lincar form of 
the Euler equations with respect to the primitive flow 
variables. This approach results in an equation similar 
to that of Eq. (9) given by 




A(?j 


m= 1 

4 


+ ^2 a (gmlAS A<7,„ = 


- ^ E [ 7/ (^ + ) + 11 (<?■) - H (<? + - <n]"A.9 

m = l 

(12) 

where the matrix D relates the time derivative of Q to 
the time derivative of q as simply 


dQ n 0q 

0t dt 


(13) 


The discretization represented by Eq. (12) is more effi- 
cient than the discretization represented by Eq. (9) be- 
cause the flux Jacobians a + and a~ are simpler mathe- 
matically and therefore faster to compute than the flux 
Jacobians and A ~ . 


4.1 Gauss-Seidel Relaxation Procedure 

Direct solution of the system of simultaneous equa- 
tions which results from application of Eq. (12) for all 
tetrahedra in the mesh, requires the inversion of a large 
matrix with large bandwidth which is computationally ex- 
pensive. Instead a Gauss-Seidel (GS) relaxation approach 
is used to solve the equations whereby the summation in- 
volving Aq m is moved to the right hand side of Eq. (12). 
The terms in this summation are then evaluated for a 


given time step using the most recently computed values 
for A <j m . The solution procedure then involves only the 
inversion of a 5x5 matrix (represented by the terms in 
square brackets on the left hand side of Eq. (12) for each 
tetrahedron in the mesh. Also, although the procedure is 
implemented for application on (randomly -ordered) un- 
structured meshes, it is not a point Gauss-Seidel proce- 
dure. The method is in fact more like line Gauss-Seidel 
since the list of tetrahedra that make up the unstructured 
mesh is rc-ordercd from upstream to downstream, and the 
solution is obtained by sweeping two times through the 
mesh as dictated by stability considerations. The first 
sweep is performed in the direction from upstream to 
downstcam and the second sweep is from downstream to 
upstream. For purely supersonic flows, the second sweep 
is unnecessary. 


4.2 Point-Jacobi Relaxation Procedure 


The inner-most do-loop of the Gauss-Seidel relax- 
ation procedure docs not vectorize due to vector recur- 
rences resulting from the evaluation of Aq m using the 
most recently computed values. Hence, a two-sweep 
Point-Jacobi (PJ) type of relaxation procedure has been 
developed that fully vectorizes. It is consequently faster 
per iteration than the GS procedure, but is expected to 
have diminished stability properties since the first sweep 
ignores the A q m term altogether, and the second sweep 
uses the values of A q from the first sweep to evaluate the 
Aq m temi in the second sweep. Thus the PJ procedure 
is represented by 


first sweep: 


T7 ; + ^fl + (<7,)A.9 


A,] 1 ' = 
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(14a) 


second sweep: 
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(14b) 


4.3 Point-Implicit Procedure 

Advantages of the Gauss-Seidel and Point-Jacobi re- 
laxation procedures are that they are numerically stable 
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for reasonably large CFL numbers, even on very fine 
meshes, and consequently they enable rapid convergence 
to steady state. For unsteady applications, they allow the 
selection of the step size based on the physical problem 
rather than on numerical stability considerations. This is 
in contrast with an explicit time integration which has 
a restrictive step size for unsteady applications which is 
more severe on finer meshes. A disadvantage of the GS 
and PJ relaxation procedures, though, is that they require 
about twice the memory of an explicit time integration, 
primarily due to having to store the backward flux Jaco- 
bian a ~ . Hence, a single-sweep Point-Implicit procedure 
was developed (represented by Eq. (14a)) that docs not 
require the calculation of the backward flux Jacobian. It 
is consequently faster than the GS or PJ relaxation proce- 
dures, but is expected to have diminished stability prop- 
erties since the A q m term is totally ignored. 


5. BOUNDARY CONDITIONS 

To impose the flow tangcncy boundary conditions 
along the surface of the vehicle, the flow variables are set 
within dummy cells that are effectively inside the geom- 
etry being considered. The velocity components within 
a dummy cell, are determined from the val- 

ues in the cell j adjacent to the surface, (w,t\w)j- This 
is accomplished by first rotating the components into a 
coordinate system that has a coordinate direction normal 
to the boundary face. The sign of the velocity compo- 
nent in this direction is changed (hence imposing no flow 
through the face) and the three velocity components are 
then rotated back into the original x,y*z coordinate sys- 
tem. After considerable algebra this yields 
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the unit vector that is normal to the boundary face. Also, 
pressure and density within the dummy cells are set equal 
to the values in the cell adjacent to the surface. 

After application of the upwind -biased interpolation 
formula to determine q~ and at each face, the velocity 
components are corrected to give a "strong” implementa- 
tion of the surface boundary condition according to 


Ucorr'cted = tt - + vn y + Kin.) 

^corrected = V- n y (UTl z + V)l v + «'!!.) (16) 

WcoTT'rted = «’ - n z (un T + vn y + urn , ) 

In the far field a characteristic analysis based on Riemann 
invariants is used to determine the values of the flow 


variables on the outer boundary of the grid. This analysis 
correctly accounts for wave propagation in the far field 
which is important for rapid convergence to steady-state. 


6. RESULTS AND DISCUSSION 

To assess the efficiency of the implicit upwind solu- 
tion algorithms, calculations were performed for the Boe- 
ing 747 aircraft configuration. The results were obtained 
using the unstructured mesh shown in Fig. 1. The 747 
geometry includes the fuselage, the wing, horizontal and 
vertical tails, underwing pylons, and flow-through engine 
nacelles. The unstructured mesh for the 747 contains 
101,475 tetrahedra and 19,055 nodes for the half-span 
airplane. Also there are 4,159 nodes and 8,330 triangles 
on the boundaries of the mesh which include the airplane, 
the symmetry plane, and the far field. Steady-state cal- 
culations were performed for the aircraft at a frccstream 
Mach number M ^ of 0.84 and an angle of attack a of 
2.73°. All of the results were obtained on the Cray-2 
computer (Navier) at the Numerical Aerodynamic Simu- 
lation Facility located at NASA Ames Research Center. 

Steady-state calculations were performed first using 
the implicit Gauss-Seidcl relaxation procedure. These 
implicit results were obtained using a CFL number of 
infinity and the flux Jacobians were updated only every 
twenty iterations. Such a large value of the CFL number 
was used for the GS results since the relaxation scheme 
has maximum damping and hence fastest convergence for 
very large time steps. This is in contrast with implicit 
approximate factorization schemes which have maximum 
damping for CFL numbers on the order of 10. Results 
also were obtained for comparison purposes using an 
explicit, three-stage, Runge-Kutta time-marching scheme 
(ref. 2). These explicit results were obtained using a CFL 
number of 4.0, with residual smoothing and local time- 
stepping to accelerate convergence to steady state. 

A comparison of the convergence histories between 
explicit and implicit GS solutions is shown in Fig. 2(a). 
The "error” in the solution was taken to be the Li- 
norm of the density residual. As shown in Fig. 2(a), the 
GS solution converges faster than the explicit solution. 
The GS solution, for example, required 622 iterations 
or approximately 3,420 CPU secs (less than one hour) 
to converge to engineering accuracy, which is taken to 
be a four order-of-magnitude reduction in the solution 
error. In contrast, the explicit solution required 1,552 
iterations or approximately 11,560 CPU secs to achieve 
the same convergence. The GS relaxation procedure is 
not only faster on a per iteration basis, but provides faster 
convergence to steady state in terms of CPU time. 
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Figure 1 Surface mesh of triangles for the Boeing 747 aircraft. 
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(a) Gauss-Seidel relaxation procedure. (b) Point-Jacobi relaxation procedure. (c) Point-Implicit procedure. 

Figure 2 Comparison of implicit and explicit convergence histories for the Boeing 747 aircraft at M ^ = 0 84 and a = 2.73°. 



Figure 3 Steady pressure coefficient contours on the Boeing 747 aircraft at = 0.84 and o = 2 73°. 
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Calculations were performed next with the implicit 
Point- Jacobi relaxation procedure for the 747 aircraft. 
These calculations were unstable for CFL numbers of in- 
finity and 100 due to updating the Aq m term in a point 
Jacobi fashion rather than in a Gauss-Seidel fashion. Sta- 
ble results were obtained using a CFL number of 10 and 
local time stepping was used to accelerate convergence 
to steady state. As shown in Fig. 2(b), the PJ solution 
converges faster than the explicit solution, but not as fast 
as the GS solution convergence of Fig. 2(a). The PJ re- 
sult required 1,942 iterations or approximately 5,080 CPU 
secs to converge the solution four orders of magnitude. 
Hence, the PJ relaxation procedure, although not faster 
than either the explicit or GS procedures on a per iter- 
ation basis, provides faster convergence to steady stale 
in terms of CPU time in comparison with the explicit 
Runge-Kutta procedure. 

Calculations were performed also with the Point- 
Implicit procedure for the 747 aircraft. These calculations 
were unstable for CFL numbers of infinity, 100, and 10 
due to the approximations that are made to solve the 
implicit equations. Stable results were obtained using a 
CFL number of unity and local time stepping was used 
to accelerate convergence to steady state. As shown 
in Fig. 2(c), the PI solution converges slower than the 
explicit solution, and is also slower than the other two 
implicit procedures. Although the PI algorithm converges 
the slowest of all of the methods used in the present study, 
it is the fastest algorithm on a per iteration basis. The 
PI algorithm is about twice as fast as the GS relaxation 
procedure and is nearly three times as fast as the explicit 
Runge-Kutta method. Efforts to improve the rate of 
convergence of the PI procedure by slowly increasing 
the CFL number over the course of the calculation were 
unsuccessful. 

Finally, steady pressure coefficient contours on the 
surface of the 747 aircraft arc shown in Fig. 3. These 
results were obtained using the GS relaxation procedure 
with a CFL number of infinity. The contours indicate that 
there is a significant amount of flow compression on the 
nose of the aircraft along the inboard leading edge of the 
wing, and inside the cowl of the engine nacelles. There is 
flow expansion on the forward fuselage, on the horizontal 
and vertical tail surfaces, and on the upper surface of the 
wing terminated by a shock wave. 

7. CONCLUDING REMARKS 

The development of implicit upwind algorithms for 
the solution of the three-dimensional time-dependent Eu- 
ler equations on unstructured tetrahedral meshes was de- 
scribed. The implicit temporal discretization involves ei- 


ther a two-sweep Gauss-Seidel relaxation procedure, a 
two-sweep Point-Jacobi relaxation procedure, or a single- 
sweep Point-Implicit procedure; whereas the upwind spa- 
tial discretization is based on the flux-difference splitting 
(FDS) of Roe. The FDS discretization is naturally dissi- 
pative, and consequently captures shock waves and con- 
tact discontinuities sharply. Detailed descriptions of the 
three implicit solution algorithms were given, and calcu- 
lations for the Boeing 747 transport configuration were 
presented to demonstrate the algorithms. The 747 ge- 
ometry included the fuselage, wing, horizontal and verti- 
cal tails, underwing pylons, and flow-through engine na- 
celles. Advantages and disadvantages of the implicit al- 
gorithms were discussed. A steady-stale solution for the 
747 configuration, obtained at transonic flow conditions 
using a mesh of over 100,000 cells, required less than one 
hour of CPU time on a Cray-2 computer, thus demonstrat- 
ing the speed and robustness of the general capability. 
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